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ABSTRACT 


This paper presents error analyses for numerically 
integrating first order ordinary differential equations. 
Roundoff, truncation, mathematical instability, and numer- 
ical instability errors are the four basic errors considered. 
Emphasis is placed upon the catastrophic effects of numeri- 
cal instability errors. Several examples of numerical and 
mathematical instability errors are contained in the paper. 
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AN ENGINEERING EXPOSE OF NUMERICAL INTEGRATION 

OF ORDINARY DIFFERENTIAL EQUATIONS 

By John L. Engvall 
Manned Spacecraft Center 

SUMMARY 


This paper presents error analyses for numerically integrating first order 
ordinary differential equations. Roundoff, truncation, mathematical instability, and 
numerical instability errors are the four basic errors considered with emphasis 
placed upon the catastrophic effects of numerical instability errors. Several examples 
of numerical and mathematical instability errors are presented. 

Although instability is not a totally undesirable characteristic, in many cases 
the effect of instability errors is more significant than roundoff errors and trunca- 
tion errors. 


INTRODUCTION 


Mathematical analysis concerning numerical integration of ordinary differential 
equations has been in progress at the Manned Spacecraft Center (MSC) for some time. 
This paper presents the theoretical results of a portion of this error analysis, and 
contains numerical examples which indicate the order of magnitude of the difficulties 
that may be encountered due to various errors. Emphasis is placed upon the impor- 
tance of error growth due to numerical instability of multistep methods. 


SYMBOLS 


A 

A(x) 

C l» C 2’ C 3’ C 4 


array of coefficients for infinite series associated with numerical 
methods 

particular solution to differential equation 
constants of integration 


dy/dx = f(x, y) denotes ordinary differential equation 



E, ER(x) 

ET(x) 

exp(x) 

F(x) 

G(x) 

HBAR 

h 

P 

q 

R l’ R 2’ R 3’ R 4 
RER(x) 

RET 

T 


error in numerical solution 

local truncation error 

exponential function of x 

total solution to differential equation 

homogeneous solution to differential equation 

product of step size and 3f/ Sy 

denotes the integration step size 

present set of data 

array of differences T(I) - A(I) 

expected relative errors for example 1 

relative total error 

relative local truncation error 

array of Taylor series constant coefficients 

initial independent variable 


x^j approximate initial condition 

yN, FN numerical solution to differential equation 


y (n) « 

y (‘o) 

y(x) 

3f/3y 


denotes the nth derivative of y with respect to x 
exact initial condition for solution to differential equation 

denotes solution in differential equation above 
partial derivative of f with respect to y 


Symbols not listed are used as constant parameters in the development of 
equations. 
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BASIC DEFINITIONS 


Definition 1 

Roundoff error is any error which is directly caused by the fact that, in com- 
puters, numbers must be expressed by a finite number of significant figures. 

The term "truncation error" is often applied to error introduced by truncating a 
number or by retaining only the most significant part of a number. According to the 
above definition, this is roundoff error. The term "truncation error" will be used 
for a different type of error and is defined subsequently. 

Roundoff error can occur before what is generally considered to be initial execu- 
tion of the numerical integration. Conversion of input data from decimal to binary 
might generate errors of this form. For example, a number which can be represented 
by a finite number of significant figures in the decimal system may have an infinite 
repeating binary representation; for example, 

0. 1 (base ten) =0. 000110011001100 (base two) 

Normally, individual roundoff errors tend to occur in the least significant fig- 
ures carried within the computer. This certainly seems insignificant, and, as a 
matter of fact, roundoff error often receives little treatment in applied texts on 
numerical methods. Actually, roundoff error is the most significant error in the 
sense that if roundoff error can be made arbitrarily small, then almost all other 
errors can be made arbitrarily small; yet paradoxically, roundoff error can usually 
be assumed to be negligible compared to other errors. The above statement will be 
discussed later in detail. 

Discussion of any other type errors will be governed by the following five 
restrictions: 

a. Only differential equations which have unique bounded solutions through the 
initial data point and over the domain of numerical integration will be considered. The 
necessary and sufficient conditions for this restriction are fully discussed in refer- 
ence 1. 


b. The dependent variable is assumed to have derivatives of all orders, with 
respect to the independent variable, over the domain of any one integration step. 

c. The numerical approximation given by any method considered herein can 
theoretically be expressed as an infinite series. This infinite series differs from 
Taylor's series about a given data point only in the values of the constant coefficients. 

d. The numerical methods will be restricted to two basic types. These types 
are the single-step methods such as Runge-Kutta, and the multistep methods con- 
sisting of explicit or implicit methods utilizing linear difference equations, such as 
Adams-Moulton, using equally spaced points. Development of these methods can be 
found in reference 2.- 
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e. Only first order systems will be considered. Any errors, which might arise 
from integrating through or near a point where a unique bounded solution does not exist, 
will not be considered. The Taylor series, in restriction c above, is expanded about 
the single point for single step methods and one of the points contained in the difference 
equation for multistep methods. The series approximates the value of the dependent 
variable at a distance h away from this point. 


Definition 2 

The distance h shall be referred to as the integration step size. 


Definition 3 

A present set of data is a set of data which has been obtained by conversion from 
decimal to binary and/or integration of the differential equation by some numerical 
method from some initial set of data. 

Roundoff error usually precludes the possibility of zero error in a present set of 
data. The error can be assumed to be zero and then the error of one integration step 
with respect to a present set of data can be measured. 

Generally, the numerical methods used to solve differential equations are based 
upon a mathematical model which yields an erroneous solution for each integration 
step. One of the primary reasons for this is that many methods are based on an exact 
mathematical solution expressed by using infinite difference tables or infinite series. 
The numerical method is then based upon a finite process which differs from the true 
solution but agrees with the infinite process up to some desired term. This difference 
is usually referred to as truncation error. For the restrictive types of methods dis- 
cussed in this paper, truncation error can be described as follows. 

Given a method, a differential equation, a present set of data, and an integration 
step size, all satisfying restrictions a through e of Definition 1; then, (1) there 
exists a Taylor series expansion about the present value of the independent variable, 
and (2) there exists an expansion about the same point which is associated with the 
numerical method and differs from (1) only in the values of the constant coefficients. 


Definition 4 

Assuming that the present set of data is exact (zero error), then the local trun- 
cation error is the difference between the Taylor series expansion and the expansion 
associated with the method being used. This is a measure of the error introduced in 
a single integration step. 

Assuming the data to be exact simply means that one is defining error with 
respect to the solution through the present point rather than with respect to the initial 
point. Note that after several integration steps there might exist a large amount of 
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error in the present set of data with respect to the initial data. Thus the local trunca- 
tion error need not reflect the size of the total error with respect to the initial data. 

For the remainder of this paper the term "truncation error" refers to local 
truncation error unless otherwise stated. 


Definition 5 

The value of the independent variable associated with a present set of data shall 
be referred to as the point P. 

For example, if integrating from the point P to P + h, the truncation error for 
this step is usually approximated by using two (or more) approximations to the value 
of the dependent variable at P + h. The difference in these two values is then used as 
a criterion for determining the approximate truncation error for that step. If the 
error is acceptable, the dependent variable at the new point can be reinitialized and 
the integration procedure applied again. This is effectively the equivalent of expand- 
ing a new series about the point P + h. 


Definition 6 

The order of a method is the largest value of n for which the series associated 
with the method agrees with all terms in the Taylor series through the term containing 
the nth derivative. This is compatible with the definitions in reference 2. 

Consider the infinite series for yN(x + h) associated with a given method at a 
point x. 


yN(x + h) = S = y(x) + A(l)hy^(x) + A(2)h 2 y^(x) + . . . 


where A(I) are the coefficients associated with the given method. 

The Taylor series at the same point is 

y(x + h) = T = y(x) + T(l)hy (1) (x) + T(2)h 2 y (2) (x) + . . . 

where: 

T(J) = 1/J J = 1,2, . . . 
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Define 


q(I) = T(I) - A(I), I = 1, 2, . . . 


If a method is of order K then: 


q(i) = 0, 1 = 1,2, . . . , K 


If the truncation error is not acceptable, there is the option of decreasing the 
magnitude of h or using a higher order method. It might be expected that the numer - 
ical solution becomes exact in the limit as h approaches zero (for a given fixed 
method and order). This would be true only if an infinite precision computer were 
used. It might also be expected that for a fixed step size, the error would tend to 
decrease as higher and higher order methods are used. Although truncation error 
might be made arbitrarily small by decreasing h, there is again the impeding effect 
of roundoff error. For example, suppose h is decreased until the numerical solution 
agrees with the true solution to within the least significant figure which can be carried 
in computation (for one integration step). To reduce h beyond this point would 
decrease the truncation error, but it could not possibly increase the number of accu- 
rate significant figures. 

Truncation error can be decreased by using higher order methods; however, if 
multistep methods are used, numerical instability can cause the errors to increase 
markedly. (This will be demonstrated later. ) Higher order single-step methods (for 
example, Runge-Kutta methods) involve the use of parameters of large magnitudes 
and the calculation of more intermediate quantities. This increases roundoff error , 
and computation time. Thus, a low-order, single-step method can possibly achieve 
better accuracy than a high-order, single-step method by using a smaller step size. 
Also, the low-order method might use less computation time because of the difference 
in the number of intermediate calculations for each integration step. 


Definition 7 

The relative truncation error is the numerical value of the local truncation 
error divided by the true solution with respect to the present set of data. 

Relative truncation error is a measure of the number of significant accurate 
figures. The transition from absolute error to relative error is essential for clarity 
in some of the subsequent numerical examples. Consider a method (using 20 signifi- 
cant decimal digits) of calculating exp(x) such that the answer is always accurate 
through 19 significant figures. As x increases without bound, the absolute error 
becomes unbounded. From this viewpoint, the process is unstable. The relative 

- 18 

error is bounded in magnitude by 10 . In this sense, the process is stable. 
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At this point, the behavior of the errors discussed could be encountered in 
evaluating a definite Riemann integral. Beyond this point, the errors discussed pre- 
sent no problem in dealing with derivatives which are functions of the independent 
variable alone. 


CONCEPT OF STABILITY 


Let y(t) be the exact solution to the differential equation y' = f(y, t) with the 
initial condition y^t^ = y^ Let x(t) be an approximate solution to the same differ- 
ential equation using the approximate initial condition x ^q) = x q- 


Definition 8 

The approximate solution x(t) shall be defined as stable provided that for every 
positive number e there exists a positive number 6 such that if 


l x ( t o)-y('o)| < 6 

then 

I x(t) - y (t) I < e 


for all t> tg. 


The above definition is one of the classical definitions of stability of a solution to 
an initial value problem. If the approximate solution is stable, then given a tolerance 
e, regardless of how small, there exists a positive 6 such that if the initial condition 
x(tp) can be determined to within a distance 6 of the exact initial condition y^t^, then 

there is assurance that the error in x(t) is less than e in magnitude for all time t 
greater than t^. The following is the negation of the definition of stability. 
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Definition 9 


The approximate solution x(t) is said to be unstable provided there exists a 
positive Cq such that for every positive 6, there exists some t^ > t Q and some 

x^q) such that 

I x (*o) - K‘o)l < 6 

and 

I x (*i) - y(‘i)l ^ e o 

Stability of a solution is a desirable characteristic; however, instability is not 
necessarily an undesir able characteristic. For example, suppose |x(t) - y(t)| can be 
made arbitrarily small for all t such that tg < t < tg. The solution x(t) would be 

acceptable although it might be unstable by the above definition. It should be noted 
that some authors might define stability with respect to a finite interval rather than 
for all t > tp. 


In reality, almost all numerical solutions are unstable in the strict sense of the 
definition. For example, in computational work, all quantities are usually restricted 
to a fixed maximum number of significant figures. Thus, if y(t) = 1. 1 at t = 1, 
and a computer having a maximum word length of 27 significant binary bits is used, 

- 2 8 

then obviously, |x(l) - y (1) | > 2” if the answer x(l) is stored in a single word. 
Although this is an unstable solution, in many cases, if . five or six significant decimal 
figures can be retained then the solution is considered satisfactory. Thus, even 
though a solution is unstable, it might be satisfactory in reality. From the engineer- 
ing viewpoint, a solution may be considered stable if for a fixed tolerance Cq there is 

a particular such that | x(t) - y(t) | < e Q for every time t such that t Q < t < tj 

provided | x 

sonable interpretation of the physical problem. 

Many authors define stability in different ways; some leave the definition to the 
intuition of the reader. The term stability shall be used in this paper as formally 
defined at the beginning of this section. 

In dealing with numerical solutions to ordinary differential equations, there is 
a natural division of stability problems into two basic categories. The first deals 
with the properties of the differential equation and will be referred to as ’’dynamic 
stability” or ’’mathematical stability. ” This form of stability can be defined in terms 
of the differential equation with no reference to a numerical method. The second 
category is dependent upon the numerical method which is used to obtain the solution 
of the differential equation. It has been pointed out that, in the strict sense of the 


(*o) ~ y^o)! ^ 6 0* * n this case ’ e 0 and 6 0 cou ^ k e derived from a rea- 
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definition, any numerical solution is unstable. This is not the form of numerical sta- 
bility to be used here. Numerical stability shall be defined according to the formal 
definition of stability where the approximate solution x(t) is the solution obtained by 
a numerical method, using infinite precision numerical computation. The last part of 
this paper discusses numerical stability of a particular class of methods for solving 
differential equations. The following sections of this paper give specific examples 
which show that instability errors can become very significant over finite domains. 


A Typical Problem 

Consider the difference in two methods of approach to the following problem. 


dy/dx = f(x, y) 


y = y Q at x = x Q 


Find a function F(x) such that y = F(x) satisfies the above conditions. 

Method 1 .- Solve the differential equation for the family of solutions. 

y(x) = G(x) + CA(x) (1) 

Solve for C using the initial conditions 

c o ■ [y 0 - G ( x o)]/ A ( x o) <2) 


Thus, 


F(x) = G(x) + C 0 A(x) 


( 3 ) 


Method 2 . - Suppose a numerical integration method is used for one step of size 
h. Let us assume that after this integration step there exists a truncation error, a 
roundoff error, or both. This error, or the combination of the two errors, shall be 
designated by E. Thus the present set of data is: 

Y « f(x q + h'j + E at x = x Q + h (4) 
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For simplicity, assume that no more errors are made in the numerical integra- 
tion. The numerical solution will be equivalent to using the present set of data as 
initial conditions and solving for C in equation (1) above. Suppose C = Cl. 

*If the numerical solution is designated by FN, then 

FN(x) = G(x) + Cl A(x) (5) 

Definition 10 

ER(x) is the difference between the true solution (using the initial conditions) 
and the numerical solution of a differential equation. In this example we have 

ER(x) = F(x) - FN(x) = (Cl - C Q ) A x = CON A(x) (6) 


The major observation at this point is that only one error introduced upon the 
first integration step (or any other step) might very definitely affect the error in fol- 
lowing steps. Interest shall be in two basic types of functions A(x) : 

a. Those which tend to amplify error. 

b. Those which tend to damp out the error. 

(A single function may be of type a over one domain, and type b over another domain. ) 


Definition 11 

Error amplification of type a in Definition 10 shall be referred to as error due to 
the "inherent characteristics" of the differential equations. 


Definition 12 

Differential equations of type b in Definition 10, which are characterized by 
A(x), shall be referred to as being mathematically stable provided that A(x) — 0 as 
x becomes unbounded, and A(x) bounded for all x > Xq. 


1 FN(x) 

conditions. 


is sometimes referred to as a solution corresponding to modified initial 
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Definition 13 


RER(x) is the value of ER(x) divided by the true solution. 

RER(x) shall be referred to as the relative error (of the numerical solution) 
with respect to the initial conditions. 

EXAMPLE: 

The following examples indicate the extent to which A(x) can affect accuracy. 
Consider the differential equations: 

a. dy/dx = e 

b. dy/dx = y 

c. dy/dx = 2e x - y 

d. dy/dx = -14e x + 15y 

Consider each equation as having the same initial condition: 

y = 1 at x = 0 

The solution to each equation is the same using this initial condition. The solu- 

X 

tion is y = e . 

The general solutions are respectively: 

a. y = e + Cj 

b. y = C 2 e x 


11 





II 


I IIIIIIIIIBII 


II I I 


I ■■■ HIM ■ ■ 


II II I I 


X ^ -X 

c. y = e + C 3 e 

d. y = e x + C 4 e 15x 


with the given initial conditions: 



c 4 = 0, c 2 = 1 


x 


The differential equations were integrated numerically in such a manner that the 
relative local truncation error magnitude was less than 10 - *® for each integration 

n 

step. The equations were integrated separately using one method. The absolute 
errors in these numerical solutions are given in table I. Note that the first three 
solutions were well behaved. RER(x) is less than or equal to the allowable relative 
local truncation error per step. (See table n. ) The fourth equation introduced 
extreme mathematical instability, forcing even the relative error to grow exponen- 
tially, as would be expected. 

The relative errors which would be expected can be written as follows: 


a. R 1 = C 1 /e X 

b. R 2 = C 2 - l/e x 

c * R 3 = C 3 e 
d. R 4 = C 4 e 14x 




* 


9 -f-Tl 

E «: n shall be used to denote 10 in most of the tables. 
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EXAMPLE (MORE GENERAL) : 

Consider the differential equation 


dy/dx = (1 + B)e X - By (7) 

y = 1 at x = 0 
B = 0 


The family of solutions to the above differential equation is 


y = 


x „ -Bx 
e + Ce 


( 8 ) 


The solution with the given initial conditions is 


y = 


X 


e 


(9) 


If integrated in the positive x direction, then truncation and roundoff errors 

are damped according to CON e , where CON can be determined as in equa- 
tion (6). Also, the larger the value of B, the faster the error is damped. Thus, as 
x increases, RER(x) approaches the allowable relative local truncation error. At 
least this is true if the errors discussed so far are the only errors present in the 
numerical integration. 


Definition 14 

In integrating from the point to x. + h, ET^x.. + jJ shall designate the local 
truncation error for this step. 

If a method is of order K, then 

OO 

ET ( x i + i)= (I 

n=K+l 
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Definition 15 


The numerical solution obtained for the dependent variable at x = x^ shall be 
designated by yN (x.). 

Assume that h is held constant. A bound will be derived for RER(x) of equa- 
tion (7). 


y(x) = e x + Ce” Bx 

(ID 

y n (x) = e x + C (-l) n B n e" Bx 

(12) 

OO 


ET ( x i+ l) = S a n hny "( x i) 

(13) 


n=K+l 


OO 0O 



n=K+l n=K+l 


" Bx i 




(14) 


where: 


C (X.) = 




Since y = e x is the solution corresponding to the exact initial conditions, it 
follows that: 


V 


et ( x i+ i) = e 1 7 v * 1 " + er ra 7 (- nB V" 


n=K+l 


n=K+l 


i 


(15) 


14 


r 


Definition 16 

RET ( X i+l) sha11 designate ET ( x i+i) divided by the true solution with respect 
to the initial conditions. 

For this case, 


RET ( x i + l) = ET ( x i + l^ /e 


)/ 6 


(16) 


ret ( x 1+ i)= 2 % h " + [ ER ( x i)/ eX ‘] 2 ( - i)nB v n 


(17) 


n=K+l 


n=K+l 


and 


RET ( x i + l) = 2 V'" + RER ( X i) 2 ( - l)V V n 


(18) 


n=K+l 


n=K+l 


if B is held constant and assuming h is small enough such that the series are 
convergent, then 


^ (- 1) n B n q n h n = constant s CON 

n=K+l 

Also, 


4 


2 


n=K+l 


q^h = constant = Rq 




(Note that R Q = RER(x) for B = 0) 
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Thus, 


RET(x. +1 ) = R Q + CON RER (x.) 


(19) 


If in general, we have the true solution, yT(x), with exact initial conditions 
given by 


yT(x) = G(x) + CA(x) 


( 20 ) 


and at a point x = x^ 


y N ( x i) - G ( x i) + c ( x i) A ( x i) 


( 21 ) 


then, 


W(x, +1 ) = [C -C(x.)] A(x 1+1 )/yT (x i+1 ) + RET(x 1+1 ) 


( 22 ) 


For the example in question, 


yT(x) = e 


(23) 


A(x) = e 
C = 0 


■Bx 


3 This implies that RER(x) and RET(x) must have the same directional sense. 
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1 


Thus, 


RER( x i+ i) = 


- C ( X i) 



+ RET (X 1+l) 


(24) 


By the triangle inequality, 



-Bx. 


i+1 


-Bx, 


< 1 


-X. - -X. 

1+1 - 1 

3 = e 


x. 

yN ( x i) " e 1 = ER ( x i) 


Thus, 


c ft) 



ER (*,) 



RER (x.) 


(27) 
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Substituting equation (27) into equation (25), 


RER ( x j + i)[ S |RER (Xl) | ♦ |RET(x 1+l) | 


(28) 


Substituting from equation (11), 


RER (x. +1 ) | ^ | RER (x.) | + | RER (x.) CON + R Q 


(29) 


The final result for the error bound is 


RER ( x i + l) I S (‘ - l CON |) l RER ( x i)l + l R 0 


(30) 


From equation (19), the necessary equation to evaluate CON is available. 

CON = 


CON = 


RET(Xj) - R q]/^ ER (x Q ) 

ET ( x i)/ eXl - R o]A R ( x o) 


(31) 

(32) 


EXAMPLE: 

Suppose an error was purposely introduced into the initial conditions such that: 


RER(x q ) = RER(O) = 20 


This is a 2000 percent error. 


yN(0) = 20 = e° + C (x Q ) e° 


x 0 = 0 


x = 19 


c ( x o) 
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Thus the calculation of CON can be obtained by the following substitution into 
equation (31). 



For multistep methods, numerical instability errors are sometimes present. 
Numerical stability theory predicts exponential error growth if one integrates outside 
the region of stability. 

The preceding theory has been applied to integration with a fourth order 
Runge-Kutta method. The numerical solution behaves exactly as predicted. This does 
not mean that the relative error is small; with a large value of B, the value for CON 
is very large. However, according to equation (30) the magnitude of the relative error 
at x^ + j is bounded by a linear function of the error at x. and the theoretical results 

provide an error bound. 

Actually, the behavior of the numerical solution can be reasonably predicted 
from knowledge of the first integration error alone. The results are quite interesting 
but too lengthy to include. 

Table HI demonstrates the results of integrating the equation over 70 integration 
steps using Runge-Kutta fourth order with h = 0. 1 for each step. 

The third column lists the results of integrating the equations. 


B = 1 y f = 2e x - y 

B = 2 y* = 3e X - 2y 

B = 19 y» = 20e X - 19y 


Initial conditions for each equation were 


y = 1 at x = 0 
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I I I 


II I III ■ 


I III 


I nil 


The fourth column lists the results of integrating the same equations but using 
erroneous data for initial conditions. 


y = 20 at x = 0 thus RER(O) = 20 


Note that the errors were completely damped. 


STABILITY OF MULTISTEP METHODS 


Recent analysis has shown that multistep methods might be very undependable 
for integrating a system of simultaneous equations. 

Numerical stability has more than one definition and many modes of application. 
The emphasis is to be placed upon the catastrophic errors rather than their cause. 
However, we shall point out that for multistep methods there is a direct correspond- 
ence between stability of the numerical solution and HBAR. The error of the numeri- 
cal solution of the point t = t Q + mh has a term E(t) = cA m where c is a constant, 

and A is a function of HBAR. The "stability of the method" is defined using this 
error term if | A | is, respectively, less than, equal to, or greater than 1, then this 
error term tends to, respectively, zero, remains constant, or becomes unbounded 
as t increases without bound. The method is called strongly stable, weakly stable, 
or unstable. It should be pointed out that this does not imply that the numerical so- 
lution is stable. The total error is also determined by truncation error, mathematical 
stability errors, and perhaps many other factors all of which are not included in the 

term cA m . There are numerous texts and journal publications which give mathemati- 
cal derivations for the case of a single differential equation (refs. 3, 4, 5, and 6). 

All of the present theory is based on linear error theory; thus, the example will be a 
linear equation which agrees exactly with the theory. 


EXAMPLE: 

Consider the previous differential equation 


dy/dx = (1 + B) e x - By 
y = 1 at x = 0 


(33) 


Define 


HBAR = 



20 



so that for 


B = 1, 2, 3, . . . and h = 0. 1 


then, 


HBAR = -0. 1, -0. 2, -0. 3, . . . 


Stability theory predicts that numerical stability is a function of HBAR. The 

4 

following tables indicate the sensitivity of Adams-Moulton methods to the value of 
HBAR. Table IV shows results of integrating 70 steps with h = 0. 1. Table V shows 
the results of integrating with h = 0.01. The results of the fourth order 
Adams-Moulton and the fourth order Runge-Kutta were compared. Both results were 
obtained with the same step size and the same number of integration steps. 

Table VI is a more dynamic example of the catastrophic effect encountered when 
numerical instability errors occur. RER(7. 0) can be obtained approximately by 

3 

dividing through by 10 . 

The accuracy of the solutions in tables IV, V, and VI could have been improved 
by applying the Adams-Moulton corrector more than once at each integration step and 
successively substituting the corrected value of y into the corrector equation. The 
following equation was integrated from x = 1. 5 to x = 1. 6 with one integration step 
of h = 0. 1. 


dy/dx = lle x - lOy 


with initial condition 


y = 20 at x = 0 


The solution is: 


/ \ x 1fl -lOx 
y(x) = e + 19e 


^Adams-Moulton corrector, Adams-Bashforth predictor - both the same order. 
Initial conditions were accurate to within 15 significant decimal digits. 
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The numerically integrated value for y(l. 6) was substituted back into the cor- 
rector equation for 10 iterations. All input points were exact to within roundoff error. 
Table VII gives the results for four different orders of Adams-Moulton integration. 

Even after 10 iterations the sum of truncation and instability errors for the 
14th order and 15th order methods were greater than that of the 3rd order method. 
Note that this is for one integration step. The error would grow much faster for the 
high -order methods as the number of integration steps increased. Table Vm shows 
the result of integrating the following equation from x = 0 to x = 7. 0 with h = 0.1. 


dy/dx = 20e x - 19y 


With initial condition 


y = 1 at x = 0 


Thus, 


y(x) = e x + Ce"^ x , C = 0 


Table VIII also shows a comparison for 1, 2, and 10 iterations. 

If the corrector is used iteratively, then the difference between consecutive 
answers should be used as a criterion for determining the number of iterations. 
Choosing the number of iterations arbitrarily can result in hazardous errors as shown 
in column three of table VIII. This type of behavior was also present when four itera- 
tions were used at each integration step, but the errors using three iterations per 
step were about the same as using one iteration. 


DISCUSSION AND RESULTS 


In the absence of roundoff error, the error for any integration step could be 
made arbitrarily small. The truncation error could be made small by decreasing step 
size or by using higher order methods. If roundoff and truncation errors could be 
made arbitrarily small, then so could errors due to numerical and mathematical 
instability. 

If mathematical instability errors are present, then the truncation error must 
be kept to a minimum at each integration step. If the step size is decreased, then 
more integration steps must be used thus yielding more roundoff error. If higher 
order methods are used, the coefficients in the integration formulas become larger 
thus causing more roundoff error. Roundoff error does not cause truncation error, 
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but it prevents effective control of mathematical instability by means of controlling 
truncation error. 

If numerical instability errors are present, then a lower order method can be 
used or the step size can be decreased. If a lower order method is used, then the 
step size must be decreased to retain the same truncation error. Decreasing step 
size increases the number of integration steps which increases roundoff error. 


CONCLUDING REMARKS 


A low-order single-step method requires the use of a smaller integration step 
size than a higher order method if the truncation errors are to be the same for both 
methods. If a higher order multistep method is used, then numerical instability 
errors might cause the solution to be less accurate than that of the low-order method. 
Higher order single-step methods increase the size of the constant parameters used 
in computation and the number of times that the derivative must be calculated. This 
can increase roundoff error and the computation time required for the integration. 

If HBAR is calculated at each integration step to eliminate numerical stability 
errors, then a higher order multistep method could be used successfully. However, 
for large systems of equations, this calculation becomes time consuming as it 
involves solving for the eigenvalues of a matrix and the zeros of polynomials. (For a 
development of this theory, see reference 7. ) Thus, the use of high-order multistep 
methods introduces serious difficulty. 


Manned Spacecraft Center 

National Aeronautics and Space Administration 
Houston, Texas, August 19, 1966 
030-00-00-CF-72 
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TABLE I.- ABSOLUTE ERRORS 

Initial conditions for all equations: y = 1; x = 0 
Solution to all equations: y = e x 


Independent 

Absolute errors 

variable x 

Equation 1 

Equation 2 

Equation 3 

3.0 

0. 134 E - 10 

0.603 E - 11 

0. 538 E - 10 

6.0 

. 283 E - 09 

. 734 E - 09 

. 108 E - 08 

9.0 

. 571 E - 08 

. 271 E - 07 

. 217 E - 07 

12.0 

. 114 E - 06 

.791 E - 06 

. 435 E - 06 

15. 0 

. 230 E - 05 

. 203 E - 04 

.875 E - 05 

18.0 

.463 E - 04 

. 518 E - 03 

. 176 E - 03 

21.0 

.929 E - 03 

. 124 E - 01 

. 353 E - 02 

24.0 

. 187 E - 01 

. 289 E - 00 

. 709 E - 01 

Equation 4 

0.752 

0. 309 E - 07 



1. 5 

. 238 E - 02 



2. 03 

.621 E + 01 



2. 36 

.875 E + 03 



3. 01 

. 165 E + 08 



4. 04 

.861 E + 14 



4. 51 

.975 E + 17 



5.07 

.450 E + 21 



5. 54 

. 509 E + 27 



6.0 

. 576 E + 27 
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TABLE n. - RELATIVE ERRORS 


Independent 
variable x 

Relative errors 

Equation 1 

Equation 2 

Equation 3 

3.0 

0.659 E - 12 

0. 294 E - 12 

0. 263 E - 11 

6.0 

.689 E - 12 

. 178 E - 11 

.263 E - 11 

9.0 

.692 E - 12 

. 327 E - 11 

. 263 E - 11 

12.0 

.692 E - 12 

.477 E - 11 

. 263 E - 11 

15.0 

.692 E - 12 

.626 E - 11 

.263 E - 11 

18.0 

.692 E - 12 

. 776 E - 11 

. 263 E - 11 

21.0 

.692 E - 12 

.925 E - 11 

.263 E - 11 

24.0 

.692 E - 12 

. 107 E - 10 

. 263 E - 11 

— 

Equation 4 

0. 377 

0. 76 E - 10 



.752 

. 146 E - 07 



1.5 

. 530 E - 03 



2.03 

.818 E - 00 



2. 36 

• 8 29 E + 02 



3.01 

. 810 E + 06 



4.04 

. 151 E + 13 



4.51 

. 106 E + 16 



5.07 

. 281 E + 19 



5. 54 

. 199 E + 22 



6.0 

. 141 E + 25 



6.48 

. 999 E + 27 
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TABLE m. - MATHEMATICALLY STABLE SOLUTIONS 


[Runge-Kutta] 



B 

h • 3f/0y 

RER (7. 0) 

RER (7. 0) when introducing 
2000 percent error on 1st 
step 

1 

-0.1 

0.81 E - 06 

0. 16 E - 04 

2 

-.2 

. 35 E - 05 

. 36 E - 05 

3 

-.3 

. 88 E - 05 

. 88 E - 05 

4 

-.4 

. 17 E - 04 

. 17 E - 04 

5 

-.5 

. 28 E - 04 

. 28 E - 04 

6 

-.6 

. 43 E - 04 

. 43 E - 04 

7 

-.7 

. 62 E - 04 

. 62 E - 04 

8 

-.8 

. 86 E - 04 

. 86 E - 04 

9 

-.9 

. 11 E - 03 

. 11 E - 03 

10 

-1.0 

. 15 E - 03 

. 15 E - 03 

11 

-1. 1 

. 19 E - 03 

. 19 E - 03 

12 

-1. 2 

. 23 E - 03 

. 23 E - 03 

13 

-1. 3 

. 29 E - 03 

. 13 E - 03 

14 

-1.4 

. 36 E - 03 

. 36 E - 03 

15 

-1. 5 

.43 E - 03 

.43 E - 03 

16 

-1.6 

. 53 E - 03 

. 53 E - 03 

17 

i 

-1.7 

.63 E - 03 

.63 E - 03 

18 

-1.8 

. 76 E - 03 

. 76 E - 03 

19 

-1.9 

.92 E - 03 

. 92 E - 03 
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TABLE IV. - ADAMS-MOULTON, ADAMS- BASH FORTH 
[RER(x) evaluated at x = 7. 0 after 70 integration steps of h = 0. l] 


HBAR 

1st order 
RER (x) 

2nd order 
RER (x) 

3rd order 
RER (x) 

4th order 
RER (x) 

0.0 

0.83 E - 03 

0. 40 E - 04 

0. 24 E - 05 

0. 16 E - 06 

1 

. 54 E - 03 

. 28 E - 04 

. 18 E - 04 

. 13 E - 06 

-.2 

.45 E - 03 

. 24 E - 04 

. 16 E - 05 

. 12 E - 06 

-. 3 

. 41 E - 03 

. 23 E - 04 

. 16 E - 05 

. 12 E - 06 

-.4 

. 40 E - 03 

. 23 E - 04 

. 16 E - 05 

. 12 E - 06 

-. 5 

.40 E - 03 

. 23 E - 04 

. 16 E - 05 

. 12 E - 06 

-.6 

. 40 E - 03 

. 23 E - 04 

. 16 E - 05 

. 13 E - 06 

-.7 

.42 E - 03 

. 24 E - 04 

. 17 E - 05 

. 13 E - 06 

-.8 

. 44 E - 03 

. 25 E - 04 

. 18 E - 05 

. 14 E - 06 

-.9 

. 47 E - 03 

. 26 E - 04 

. 18 E - 05 

. 14 E - 06 

-1.0 

. 50 E - 03 

. 28 E - 04 

. 19 E - 05 

. 15 E - 06 

-1.1 

. 55 E - 03 

. 29 E - 04 

. 20 E - 05 

. 64 E - 07 

-1.2 

.60 E - 03 

. 32 E - 04 

. 22 E - 05 

. 79 E - 05 

-1. 3 

. 68 E - 03 

. 34 E - 04 

. 23 E - 05 

. 38 E - 03 

-1.4 

.78 E - 03 

. 37 E - 04 

. 24 E - 05 

. 13 E - 01 

-1. 5 

. 93 E - 03 

. 41 E - 04 

. 78 E - 05 

. 32 E - 00 

-1.6 

. 11 E - 02 

. 46 E - 04 

. 12 E - 03 

. 57 E + 01 

-1.7 

. 15 E - 02 

. 53 E - 04 

. 18 E - 02 

. 71 E + 02 

-1.8 

. 22 E - 02 

. 61 E - 04 

. 71 E - 03 

. 54 E + 03 

-1.9 

.41 E - 02 

.67 E - 04 

. 34 E - 00 

. 35 E + 03 
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TABLE V. - ADAMS-MOULTON, ADAMS- BASHFORTH 


HBAR 

12th order 
RER(. 7) 

13th order 
RER(. 7) 

14th order 
RER(. 7) 

15th order 
RER(. 7) 

0.0 

0. 11 E - 15 

0. 99 E - 15 

0. 18 E - 14 

0. 35 E - 14 

1 

• 

o 

. 12 E - 14 

.77 E - 15 

. 13 E - 14 

. 30 E - 14 

-.02 

. 12 E - 14 

. 00 E - 15 

. 99 E - 15 

. 30 E - 13 

-.03 

. 14 E - 14 

.HE- 15 

. 27 E - 14 

.99 E - 13 

-.04 

. 14 E - 14 

. 11 E - 15 

. 13 E - 14 

. 16 E - 12 

-. 05 

. 12 E - 14 

. 11 E - 15 

. 17 E - 14 

. 22 E - 15 

-.06 

.99 E - 15 

. 55 E - 15 

. 11 E - 15 

. 19 E - 14 

-.07 

. 14 E - 14 

. 44 E - 15 

. 66 E - 15 

.63 E - 13 

-.08 

. 16 E - 14 

.77 E - 15 

. 14 E - 14 

. 17 E - 11 

-.09 

. 12 E - 14 

. 38 E - 14 

. 12 E - 12 

. 14 E - 11 

1 

• 

H- 1 

O 

. 55 E - 15 

. 95 E - 14 

. 19 E - 11 

. 11 E - 09 

-. 11 

. 25 E - 14 

. 10 E - 12 

. 15 E - 10 

. 31 E - 09 

-. 12 

. 12 E - 13 

.65 E - 12 

. 32 E - 10 

. 78 E - 08 

-. 13 

. 29 E - 13 

. 20 E - 12 

. 31 E - 09 

. 38 E - 07 

-. 14 

. 26 E - 12 

. 10 E - 10 

. 27 E - 08 

. 59 E - 07 

-. 15 

. 17 E - 11 

.40 E - 10 

. 51 E - 08 

. 12 E - 05 

-. 16 

.72 E - 11 

. 49 E - 09 

. 56 E - 07 

. 64 E - 05 

-. 17 

. 30 E - 12 

. 87 E - 09 

. 50 E - 06 

. 26 E - 04 

18 

. 59 E - 10 

. 71 E - 08 

. 15 E - 05 

. 64 E - 04 

-. 19 

. 42 E - 09 

1 

. 15 E - 07 

. 43 E - 05 

. 74 E - 04 


29 


TABLE VI. - 14th ORDER ADAMS-MOULTON VERSUS 2nd ORDER 

ADAMS- MOULTON 

Both sets are using same step size, same initial conditions,! 

same number of integration steps J 


HBAR 


14th order 
ER (7. 0) 


2nd order 


ER(7. 0) 


0. 37 E - 11 
. 10 E - 10 
.41 E - 03 
. 19 E + 02 
.43 E + 05 
. 16 E + 08 
. 18 E + 10 
.40 E + 11 
. 34 E + 13 
.20 E + 15 
. 38 E + 16 
. 94 E + 16 
. 53 E + 18 
. 73 E + 19 
. 27 E + 20 
. 18 E + 21 
. 23 E + 22 
.81 E + 22 
. 10 E + 23 
. 18 E + 24 


0. 44 E - 01 
.31 E - 01 
. 27 E - 01 
.25 E - 01 
.25 E - 01 
.25 E - 01 
.25 E - 01 
.26 E - 01 
.27 E - 01 
.29 E - 01 
.30 E - 01 
. 32 E - 01 


.45 E - 01 
. 51 E - 01 


.60 E - 01 
.74 E - 01 




TABLE VH. - ADAMS-MOULTON, ADAMS- BASHFORTH INTEGRATION 
WITH ITERATIONS ON THE CORRECTOR 
[|ER(x)j at x = 1. 6 h = 0. 1 one step J 


Number of 
iterations 

2nd 

order 

3rd 

order 

14th 

order 

15th 

order 

1 

0.87 10" 4 

0. 37 10“ 5 

0.50 10“ 3 

0.82 10“ 3 

2 

. 17 10" 4 

.45 10“ 6 

. 11 10“ 3 

. 19 10“ 3 

3 

. 26 10" 4 

. 10 10“ 5 

.41 10“ 4 

.67 10“ 4 

4 

.81 10" 5 

.51 10“ 6 

.41 10" 7 

.45 10" 6 

5 

. 15 10“ 4 

.73 10“ 6 

. 11 10" 4 

. 17 10“ 4 

6 

. 13 10" 4 

. 64 10" 6 

.79 10“ 5 

. 12 10“ 4 

7 

. 14 10“ 4 

.67 10“ 6 

.87 10“ 5 

. 13 10“ 4 

8 

. 13 10“ 4 

.67 10" 6 

.85 10“ 5 

. 13 10" 4 

9 

. 14 10“ 4 

.67 10“ 6 

.85 10“ 5 

. 13 10“ 4 

10 

. 13 10“ 4 

.67 10“ 6 

.85 10“ 5 

. 13 10“ 4 
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TABLE vm. - ADAMS- MOULTON, ADAMS-BASHFORTH INTEGRATION 
WITH ITERATIONS ON THE CORRECTOR 
Jr.ER(x) at x = 7.0-h = 0. 1- 70 stepsj 


Order 

One 

iteration 

Two 

iterations 

Ten 

iterations 

1 

0. 42 10" 2 

0. 84 10 15 

0. 21 10" 2 

2 

. 67 10" 4 

. 16 10 19 

00 

O 

h- ^ 

o 

1 

4^ 

3 

. 34 10° 

. 31 10 22 

. 12 10" 5 

4 

.85 10 3 

. 10 10 25 

. 11 10" 2 

5 

. 21 10 7 

. 82 10 26 

. 19 10 1 

6 

. 66 10 10 

. 19 10 28 

. 12 10 3 

7 

. 19 10 12 

. 17 10 29 

. 29 10 5 

8 

14 

. 36 10 A * 

. 67 10 29 

. 90 10 6 

9 

. 34 10 15 

. 14 10 30 

. 16 10 8 

10 

. 15 10 16 

. 17 10 30 

. 99 10 8 

11 

. 35 10 16 

. 60 10 29 

. 16 10 9 

12 

. 53 10 17 

. 22 10 30 


13 

. 18 10 18 

. 36 10 30 


14 

. 18 10 20 

. 30 10 31 


15 

. 31 10 20 

. 79 10 32 
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"The aeronautical and space activities of the United States shall he 
conducted so as to cofjtribuie ... to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof ” 

— National Aeronautics and Space Act of 1958 
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